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Abstract 

We consider an initial-boundary value problem for a generalized 2D time-dependent Schrodinger 
equation on a semi-infinite strip. For the Crank-Nicolson finite- difference scheme with approximate 
or discrete transparent boundary conditions (TBCs), the Strang-type splitting with respect to the 
potential is applied. For the resulting method, the uniform in time L^-stability is proved. Due to 
the splitting, an effective direct algorithm using FFT is developed to implement the method with the 
discrete TBC for general potential. Numerical results on the tunnel effect for rectangular barriers are 
included together with the related practical error analysis. 
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1 Introduction 

The time-dependent Schrodinger equation with several variables is important in quan- 
tum mechanics, atomic and nuclear physics, wave physics, nanotechnologies, etc. Often 
it should be solved in unbounded space domains. In particular, the generalized 2D time- 
dependent Schrodinger equation with variable coefficients on a semi-infinite strip appears 
in microscopic description of low-energy nuclear fission dynamics j4[ [13] . 

Several approaches are developed and studied to solve problems of such kind, in par- 
ticular, see [iiai3l[6llllia[17lll8l|20l[25^ One of them exploits the so-called discrete 
transparent boundary conditions (TBCs) at artificial boundaries Its advantages are 
the complete absence of spurious refiections, reliable computational stability, clear math- 
ematical background and rigorous stability theory. 

The Crank-Nicolson finite-difference scheme with the discrete TBCs in the case of a 
strip or semi-infinite strip was studied in detail in [310I8]. But the scheme is implicit so 
that solving of a specific complex system of linear algebraic equations is required at each 
time level. Efficient methods to solve such systems are well developed by the moment in 
real but not complex situation. 

The splitting technique is widely used to simplify solving of the time-dependent Schro- 
dinger and related equations, in particular, see O Hi] [12l [HI [15l [1^. We apply the 
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Strang-type splitting with respect to the potential to the Crank-Nicolson scheme with 
rather general approximate TBC in the form of the Dirichlet-to-Neumann map. For the 
resulting method, we prove the uniform in time L^-stability under a condition on an 
operator S in the approximate TBC. 

To construct the discrete TBC, we consider the splitting scheme on an infinite mesh in 
the semi-infinite strip. Its uniform in time L^-stability together with the mass conservation 
law are proved. We find that an operator 5ref in the discrete TBC is the same as for the 
original Crank-Nicolson scheme in [7J, and it satisfies above mentioned condition so that 
the uniform in time L^-stability of the resulting method is guaranteed. The operator 5ref 
is written in terms of the discrete convolution in time and the discrete Fourier expansion 
in direction y perpendicular to the strip. Due to the splitting, an effective direct algorithm 
using FFT in y is developed to implement the method with the discrete TBC for general 
potential (while other coefficients are ^-independent). 

The corresponding numerical results on the tunnel effect for rectangular barriers are 
presented together with the practical error analysis in C and norms confirming the 
good error properties of the splitting scheme. 

Notice that the results are rather easily generalized to the case of a multidimensional 
parallelepiped infinite or semi-infinite in one of the directions. 

2 The Schrodinger equation on a semi-infinite strip and the 
spHtting in potential Crank-Nicolson scheme with an approx- 
imate TBC 

Let us consider the generalized 2D time-dependent Schrodinger equation 

ihpDt^l^ = {Ho + V)^l; for (x, y) e t > 0, (2.1) 
where Q := (0, oo) x (0, Y) is a semi-infinite strip, involving the 2D Hamiltonian operator 

:= -y [D,{BnD,^l;) + D,{BuDyij) + Dy{B2rD,^) + Dy{B22Dy^)] , 

with a symmetric and positive definite uniformly in Vt matrix of the real coefficients 
{Bpq{x^y)}'^^^^^^ and the real coefficients p{x^y) and V{x^y) (the potential) such that 
p{x^y) ^ p > in f]. Also i is the imaginary unit, > is a physical constant, Dt^ 
and Dy are partial derivatives, and the unknown wave function = ip{x^t) is complex- 
valued. 

Impose the following boundary condition, the condition at infinity and the initial one 
'0(•,^)|a^7 = 0, ||'0(x, •,t)||L2(o,y) ^ as x ^ +oc, for any t > 0, (2.2) 

ij\t=o = i^'{x,y) in (2.3) 

We also assume that 

Bn{x, y) = Bioo > 0, B^ix, y) = B2i{x, y) = 0, B22{x, y) = B200 > 0, 

p(x,2/)=poo>0, y(x,2/) = Kx), ^\x,y) = on ^X^Xo. (2.4) 



2 



for some Xq > 0, where Qx '= (0, X) x (0, Y). It is well-known that solution to problem 
( 2.1[ )-( [Z4 ) satisfies a non-local integro-differential TBC for any x = X ^ Xq (for example 
see [7J) which we do not reproduce here. 

Introduce a non-uniform mesh uJ^ in t on [0, oc) with nodes = to < • • • < < • • 



^ oc as m ^ oc, and steps := — tm-i- Let t, 



m-i/2 = ^^^^-^ and := uJ^\ {0}. 
In the differential case, the splitting in potential method can be represented as follows: 
three problems are solved sequentially step by step in time 



ihpDtijj = {AV)ijj on Q X {tm-i,t 



t — tn 



\t=tm-l 1 

ihpDtij = (Ho + V)^ on nx ij\t=tm-i = ^\t=tm-i/2^ 

^|aQ = 0, \\^jj{x,-,t)\\L2^ox) as x oo, for t e {t^-i,t^]] 



ihpDtip = {/W)tp on Q X (t 



m— 1/25 ^m] 5 



^1 



t=t. 



m-1/2 



t — tn 



for any m ^ 1, where := V — V and V{x) is an auxiliary potential satisfying 

V{x) = Voo on [Xo,oo). 



(2.5) 
(2.6) 
(2.7) 
(2.8) 
(2.9) 

(2.10) 



In the simplest case, V{x) = Voo- But, in particular, to generalize results to the case of 



a strip and different constant va! 
non-constant V] see also Section 



ues V±oQ of V{x^y) at x ^ ±(X), it is necessary to take 
below. 



The Cauchy problems (2.5) and ( |2.8 ) can be easily solved explicitly, in particular 



r, 



—I 



t=t„ 



2hp ^1 



t=t„ 



(2.11) 



Equation in (2.6) is the original equation (2.1) simplified by substituting V for V. Also 
lip and are auxiliary functions and ip is the main unknown one. This is a version of the 
Strang-type splitting [15] (though the original Strang splitting [19] was suggested with 
respect to space derivatives for the 2D transport equation; note that sharp error bounds 
for the Strang splitting for the 2D heat equation can be found in [23j). The symmetrized 
three-step form of this splitting ensures its second order of approximation for ip\t=tm with 
respect to r^. 

We turn to the fully discrete case. Fix some X > Xq and introduce a non-uniform mesh 
oJh.oo in X on [0, oo) with nodes = Xo<---<xj = X<... and steps hj := Xj — 
such that xj_2 ^ Xq and hj = h = hj for j ^ J. Let (jOh.oo '= ^h,oo\ {0}, uJh •= {xj}j=o 
and ujh := uJh\{0,X}. 

We define the backward, modified forward and central difference quotients as well as 
two mesh averaging operators in x 



Wi - W. 



h 



'j+l/2 



i-1 



2h 



■i+l/2 



3 ■ 



2h 



'i+l/2 
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where := ^^^^^+^ . 

We define two mesh counterparts of the inner product in the complex space L^(0,X): 

and the associated mesh norms || • ||^^ and || • ||^^ (of course, for mesh functions respectively 
defined on ujh or defined on uJh and equal zero at Xq = 0). Hereafter z*, Rez and Imz 
denote the complex conjugate, the real and the imaginary parts of z G C. The above 
averaging operators are related by an identity 

J 1 
(^-^' U)^. = E ^j('^^U;)h, - - {W.U^h + WjU}hj). (2.12) 

We also introduce a non-uniform mesh uJs in on [0, Y] with nodes = yo < • • • < 
yx = Y and steps Sk := yk — Vk-i- Let ujs := cJ(5\{0, ¥}. We define the backward and the 
modified forward difference quotients together with two mesh averaging operators in y 

o 7-7 Uk — Uk-1 ^ jj Uk+i — Uk _ jj Uk-i + Uk ^ jj 5kUk + Sk+iUk+i 

OyUk '— ^ , OyUk '— ^ , SyUk — , SyUk '— ^ , 

Ok ^/c+1/2 ^ ^Ok+1/2 

S -\-6 ° 

where •= Let H{uJs) be the space of functions U\ uJs ^ C such that 

U\k=o,K = 0, equipped with the inner product 

K-l 

iU,WU := 5] f/,W^;4+i/2 

k = l 

and the associated norm || • ||^^. 

We define the product 2D meshes cJh,cx) •= ^h,oo xuJs on Q and uJi^ := uJh x on Qx 
as well as their interiors uJh.oo '= ^h,oo x and c^h := ^/i x ^5- Let Fh = {(0,?/^^), 1 ^ 
k ^ K — 1} U {(xj, 0), (xj, y), ^ j ^ J} be a part of the boundary of cJh- 

Let A_jk '= ^(xj-i/2, 2/^-1/2)5 for all the coefficients A = p^Bpq^V^ with Xj_i/2 '•= 
xj-i+xj yj._i/2 '= • exploit the 2D mesh Hamiltonian operator 



dx{Biii,d^W) + d:,Sy{Bi2i,s^dyW) + s^dy{B2ihdxSyW) + dy{B22hdyW) 



where the coefficients are given by formulas = 'SyBu^^ B22h = ^x^22,-, ^i2h = 

B2ih = B12-. We also set ph = ^x^yP-^ = ^x^yV- and Vh = Actually this 

finite-difference discretization is a simplification of the bilinear finite element method for 
the rectangular mesh cJh (conserving, in particular, its L^(f]) and H^{Q) optimal error 
bounds), see Some other operators T^oh could be also exploited. 

We define also the backward difference quotient and an averaging in time 

Q^^rn ._ ^ - 1 
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M 

m=l ' 



and set ujI^ := {trn} 

The following Crank-Nicolson-type scheme with general approximate TBC was studied 

mm 



ihp^dt^^ = (T^oh + Vh)si*"^ on u^, 



Irh 



0, 



(2.13) 
(2.14) 



ihp^dt^ + ( Y B2oodydy - Voo ) st"^ 



Sioo<S"^*y on ujs, 



^0 = ^0 on cJh, 



(2.15) 
(2.16) 



for any m ^ 1. Here the boundary condition (2.15) is the approximate TBC, with a linear 



operator acting in the space of functions defined on ujsXlaj^^ and = {^^j., . . . , ^y.}- 
Also ^^hjTc ^ ^^{^j^Vk) (for definiteness) and thus = 0; we assume that the 

conjunction condition ^^^Irh = is valid as well. 



Recall that the left-hand side in the approximate TBC (2.15) has the form of the well 



known 2D second order approximation to ^Bi^Dx in the Neumann boundary condition 
(exploiting an 8-point stencil in all the directions x, y and t). 

We write down the following Strang-type splitting in potential for the Crank-Nicolson 
scheme ( [2T3l )-( [2l6l ) 



= AVh on {ujh U xj) X Us, 



5^ — ^ 



= (^oh + V'ft) on ojh, 

_|_ 

AVh on {ujh^xj) X ujs, 



*"^|r, = 0, ^'"Ir, =0, *'"|rH = 0, 



(2.17) 
(2.18) 

(2.19) 
(2.20) 



^2 ^m^^m ^ 



Dr(*,*)7 



B2oodydy — VoQ 



for any m ^ 1, where AVh V/^. 



(2.21) 
(2.22) 



5 



Obviously equations (2.17) and (2.19) are reduced to the explicit expressions 



m^TfTn—l 



-AVh 



1 + i 



4^ph 



-, on [ujh^xj) xujs. (2.23) 



The main finite-difference equation (2.18) is similar to the original one (2.13) simplified by 
substituting Vh for Vh- Here ^ and ^ are auxiliary unknown functions and ^ is the main 
unknown one. We have got the approximate TBC (2.21) by substituting respectively 



Trr 



for St^, dt^ and *y in the approximate TBC (2.15); but notice that since AVh|j=j = 0, 
actually 

_ ^m-l^ _ m ^ 1. (2.24) 

Clearly the constructed splitting in potential scheme can be also considered as the 
Crank-Nicolson-type discretization in time and the same approximation in space for the 
above splitting in potential differential problem (2. 5)- (2. 11). 



Also note that inserting formulas ( |2.23[ ) into equation (2.18) and excluding the auxiliary 



functions leads to the following equation for ^ 



ihph 



or, in another form 
r, 



i^Ph - Y(%h + Vh)\ (f'")-^*"^ = [ihpi, + -fiUoi, 



Vh) 



m<Tfm—l 



(2.25) 



(2.26) 



Equation (|2.25D can be considered as a non-standard discretization for the Schrodinger 
equation ( |2.1[ ) whereas equation (2.26) can be viewed as a specific symmetric approximate 
factorization [21j with respect to the potential of the Crank-Nicolson equation (2.13). 

To study stability for the splitting finite-difference scheme in more detail, we replace 
(2.18) and ( |2.21[ ) by their generalized versions 



i/iph = ('Hoh + Vh) ^ h F on cjh, 



on uj,, 



(2.27) 
(2.28) 



for Sill m ^ 1, where the perturbation F is given on cjh after setting F = on Fh- 



Since the Cauchy problem s (|2.5| ) and (|2.8|) do not need necessarily a discretization in 

A\4 



time, to cover both formulas (2.11) and (2.23), below we also admit an expression 

— t — 

gm _ g 2hph 



(2.29) 
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in (2.23). Obviously in both cases 5^ = 5* and |5| = 1. 



We introduce two mesh counterparts of the inner product in the complex space L (Qx)'- 

J-l K-l K-l 



(U, W)^^ := E E UjkW*,h,^^/2Sk+^/2, {U, W)^^ := {U, W)^^ + E Uj^W},- 4+1/2 

j=l k=l k=l 



and the associated mesh norms || • H^^^^ and || • ||^^. 
Proposition 2.1. Let the operators satisfy an inequality ^ 

M 

Im E ('5"^*", ^t*"')., ^0 foranyM^l, 



(2.30) 



m=l 



k=0,K 



0, where 



for any function ^: ujs x uj^ ^ C such that $° = and $| 
{$-^, . . . , $^}. Then for a solution of the splitting in potential scheme (2.23)^ (2.27)^ 
( |2.20 )^ ( 2.28D and (2.22) the following stability bound holds 



2 ^ 



Tm for any M > 1. (2.31) 



Proof We take the (•, •)^^-inner-product of equation (2.27) with a function W\ cjh ^ C 



such that Vl^lrh = 0. Then we sum the result by parts in x and y (using assumptions 
(2.4) and (2.10)), apply identity (2.12) and a similar identity with respect to exploit 



the approximate TBC (2.28) and obtain an identity 



a. = \d,W^ 



+ 5i2 ( S.d, Z ) d.-SyW" 



j=l k=l 



+ B2I ( d^Sy jS^dyW + B22S:, 



dy 1 )dyW* 



2 



'j<Jk 



+ (Vh ^ ^ ^ , W)_ + {F^, W)^^ - ■—B^^[S'^^^^J, Wj. 

for any m ^ 1, see [7J for more details. 



(2.32) 



The sesquilinear form on the right-hand side containing the five terms with coefficients 



Bpq and Vh is Hermit ian-symmetric. Thus choosing W 
imaginary part of the result, we get 



and separating the 



Owing to (2.20), (2.23) and (2.29) we have the pointwise equahties 



(2.33) 



Also taking into account equahties (2.24), we further derive 



-jTn II 2 

^ / m. 



Im F 



Multiplying both sides by and summing up the result over m = 1, . . . , M, we obtain 



M 



m=l 



M 



0||2 



^ImfF 



171=1 



(2.34) 



Applying inequality (2.30) and then equalities (2.33), we get 



M 



171=1 



|v^*"llL<llvS*°IL-. 



Tm[ max k/ph^^ - + ^SiX k/ph^^ 



M 



<llvS*T-. + iE 



m=l 



Tm max L^^^"^ . 



This directly implies bound (2.31) 



□ 



(2.23), (2.27), (2.20), (2.28) and (2.22) is uniquely solvable. 



Corollary 2.1. Let condition (2.30) be valid. Then the splitting in potential scheme 



In particular, the splitting scheme (2. 17) -(2.22) is uniquely solvable, and its solution 
satisfies an equality 



max ||A/ph^^^ 



l^h 



(2.35) 



Proof. The unique solvability follows from a priori bound (2.31), and equality ( 2.35] ) also 
is clear from (2.31) for F = 0. □ 



Remark 2.1. We emphasize that actually both the Crank- Nikols on scheme and the split- 
ting in potential scheme can be similarly considered and studied not only for the strip 
geometry but for much more general unbounded domain Q composed of rectangles with 
sides parallel to coordinate axes and having one or more separated semi-infinite strip out- 
lets at infinity. 
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3 The splitting in potential Crank-Nicolson scheme on the infi- 
nite space mesh and the discrete TBC 

To construct the discrete TBC, it is required to consider the sphtting in potential 



Crank-Nicolson scheme on the infinite space mesh for the original problem (2.1)-(2.4) on 
the semi-infinite strip 



m iTrm— 1 



— ^ 



rj2 



on cjhoo, 



Trr 



= (%h + K) J + on uj^, 



ihph = AVh on CJh,oo, 

^^L =0 ^f'^lr =0 ^^'^Ir =0 

l^h.oo ' l^h.cxo ' IJ^h.oo ' 



(3.1) 
(3.2) 

(3.3) 

(3.4) 
(3.5) 



for any m ^ 1, where Fh^oo '= ^h,cx)\^h,cx)- The perturbation is given on cJh,cx) after 
setting i^^lrhoo = 0; it is added to the right-hand side of the main equation (3.2) once 
again to study stability in more detail. 



Obviously once again equations (3.1) and (3.3) are reduced to explicit expressions 



(2.23) which are valid now on uJh.oo- Here we continue to admit expression (2.29) in (2.23) 
as well. 

Let i/h be a Hilbert space of mesh functions W\ cJh,cx) ^ C such that W\ry^ ^ = and 
YlJLi W^jkWt^ < CO, equipped with the inner product 



oo K-1 



{U,W) 



yZ ^ UjkWjj.hj^i/25k+i/2' 



j=l k=l 



Proposition 3.1. Let G i/h for any m ^ 1. Then there exists a unique solution 

to the splitting in potential scheme (3.1) -(3. 5) such that G i/h for any m ^ 0^ and the 
following stability bound holds 



2 ^ 



Tm for any M ^ 1. (3.6) 



(3.7) 



Moreover, in the particular case F = 0, the mass conservation law holds 

Proof. Given ^^^i^^ ^ there exists a unique solution G i^h to equation i\3.2\) . 
The much more general result was established in the proof of the corresponding Proposi- 



tion 3 in [23. Since expressions (2.23) are valid now on (jOh,oo^ this implies existence of a 



unique solution to the splitting scheme (3.1)-(3.5) such that G i^h for any m ^ 0. 
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We set TiohW := HohW on uJh.oo and HohW := on Fh^oo- The operator 'Hoh is 
bounded and self- adjoint in i^oh since 



{-koi^U, W)^^ = —Y,Y1 {^^^^y [(^-^) ^-^1 + Buis.dyU) d.SyW* 

3 = 1 k=l 

+ B2l(d,SyU) S.dyW* + B22S, [(dyU) dyW*] } hjS^ 

) jk 



(3i 



for any U^W G i^h, see [7j for details. 
From equation (3.2) an identity 



W 



— ? — 



Hi, 



+ Vh 



ly) +(F-,FFK 



follows, for any m ^ 1 and VF G i^h- Acting in the spirit of the proof of Proposition |2.1[ 
i.e. choosing VF = , separating the imaginary part of the result and applying 



the pointwise equalities (2.33) that are valid now on uJh,oo^ we get 



2r. 



Multiplying both sides by ^ and summing up the result over m = 1, . . . , M, we obtain 

^ M 



M||2 



0||2 



+ -5]lm(F-, 



_|_ 



m=l 



(3.9) 



□ 



The rest of the proof in fact repeats one for Proposition |2.1[ 

Corollary 3.1. Let = and = on uJh,oo\^h^ for any m ^ 1. If the solution to 
the splitting in potential scheme (3.1) -(3. 5) is such that G i/h; for any m ^ and 
satisfies an equation 



d:r. 



3=J 



III 

J on ujs-, for any m ^ 1, 



(3.10) 



with some operator = S^^^ then the following equality holds 



M 



m=l 



|^M||2 ._ ^ ||^M||2 



+ ||*fE,/i>0 for any M ^ 1. 

j=j+i 



(3.11) 



Proof Similarly to [7], equation (3.10) is equivalent to the approximate TBC (2.22) pro- 
vided that equation (3.2) is valid for j = J with F\j=j = 0. Thus the solution to the 
splitting schem e (|3.lD -(3.5) on the infinite mesh cJh,cx) is the solution to the splitting scheme 



( |2.23D , ( |2.27D , ( |2.20D , ( |2.28D and ( |2.22D on the finite mesh cJh, too. Then equality ([3Tlj) 
is obtained by subtracting (2.34) from (3.9). □ 
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Equality (3.11) clarifies the energy sense of inequality (2.30) for 5^ 



since 



inik = W^Wk + \\W\\k^\., for any W G H,,. 

By definition, the operator 5ref that has just appeared implicitly corresponds to the 
discrete TBC. Let us describe it explicitly. 

Let uJh.jo.oo '= {^j}j^jo- Then AVh = on uJh,j-i,oc x and clearly 



on (jOh,j-i,oo X for any m ^ 1. 



Therefore, if = on cJ/^,J-l,oo x ^5 and = on uJh,j,oo x ^5 for any m ^ 1, 
the splitting in potential scheme (3.1)-(3.5) is reduced on uJh,j-i,oo x uJs to the following 
auxiliary problem 



ihpocdt'^'^ = (Hoh,oo + Ko)^^*'^ on u;h,j,oc x cj^, 
^'^Ifc^cK = on uJh,j-i,oo, 
= on cJ/i,j_i,oo X ujs, 



(3.12) 
(3.13) 
(3.14) 



for any m ^ 1. Here "Hoh.oo is the limiting 2D mesh Hamiltonian operator with constant 
coefficients 



U 



Oh,. 



W ■■= — ^ ( Bx^d^d^W + B^oodydy^^ ) on a;h,oo\ c^h- 



Moreover, equation ( |3.10[ ) takes the simphfied form 

<S"'*y on a;^, for any m ^ 1. 



(3.15) 



Problem ( |3.12 )-(3.14) and equation ( 3.15] ) are the same that correspond to the original 
Crank-Nicolson scheme (2.13)-(2.16) and (after dividing (3.12) by p^o) were studied in 
detail in [7j in order to construct explicitly the discrete TBC. We recall briefiy the answer 
from [7]. To this end, introduce the auxiliary mesh eigenvalue problem in y 

-dydyE XE on cj^, ^|^^o,k = 0' E ^0. 
We denote by {£"/, Xid}^^^ its eigenpairs such that the functions are real-valued 

o 

and form an orthonormalized basis in H{ps)] here Xis > for all /. Clearly, for any 

o 

U G HiuJs)^ the following expansion holds 

K-l 

U = :F-^U^'^ := U^^^Eu where U^^^ = {J^Uf^ := ([/, Ei)^^ for 1^1^ K-l. 

1=1 

These formulas define the direct J-" and inverse J-'~^ transforms from the collection of 
values {Uk}^Si to the collection of its Fourier coefficients and back. 

In the case of the uniform mesh uJs^ i.e. Sk = S for any 1 ^ k ^ the eigenpairs are 
represented explicitly by the well-known formulas 

2 



{Ei)k := \l- sm^^ 



O^k^K, X^ 



16 • = 



- sm— , for l^l^K -I, 



2Y 
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and the transforms T and can be effectively implemented by applying the discrete 
fast Fourier transform (FFT) with respect to sines. 

Let the mesh in time uj^ be uniform. Recall that the discrete convolution of two 
sequences R^Q\ cJ'^ ^ C is defined by 

m 

[R * QY := ^ for any m > 0. 

The operator 5ref is given by a discrete convolution in t 

1 / / \ \ 

S^^^"^ = — J--^ [^Ri * j for any m > 1 (3.16) 

also involving the above transforms J-" and J-'~^ in y. Expressions for the kernel sequences 
i?/, 1 ^ / ^ K — 1, can be found in [7j, and we do not reproduce them here (see also 
p!Ol [9j for practically more convenient recurrent relations). 



Proposition 3.2. The operator Sj-ef satisfies inequality (2.30). 



Consequently, for the solution of the splitting in potential scheme (2.23)^ (2.27)^ (2.20)^ 



(2.28) and (2.22) with the discrete TBC (i.e. with S = Sj-ef), the stability bound (2.31) 
holds. 



The proof of inequality (2.30) is contained in ^J. 

If the matrix {Bp q}^^^ ^ is diagonal and ^/-independent together with p, the splitting 
in potential scheme (2T7|-(2.22) with the discrete TBC can be effectively implemented. 
Indeed, then we can apply to the main equation (2.18) and the discrete TBC (2.21) 
with S = 5ref, take into account the boundary condition ^^^|rh = (2.21) and formula 
(3.16) and obtain a collection of disjoint problems, for 1 ^ / ^ X — 1: 



ihph 



dx ( Biihd, 



+ Vi 



^jm{l) 



j=0 



= 0, 



on ujh, (3.17) 
(3.18) 



i^Poo Vi 



(3.19) 



for any m ^ 1, where 



J ' 



m{l) 



Given , the direct algorithm for computing comprises five steps. 



1. is computed on {ujh U xj) x ujs according to (2.23). 
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2. 
3. 
4. 



K-l 



1=1 



is computed by applying for any 1 ^ j ^ J. 



is computed by solving problem (3.17)-(3.19) for any 1 ^ I ^ K — 1. 



K-l 



is computed by applying ^ for any 1 ^ j ^ J. 



5. is computed on {ujh U xj) x cj^^ according to (2.23). 

Steps 1 and 5 require 0{JK) arithmetic operations, steps 2 and 4 require 0{JK log2 K)) 
operations provided that K = 2^ with the integer and step 3 requires 0{{J + m)K) 
operations. The total amount of operations is 0((Jlog2 + m)K) for computing the 
solution on the time level m and 0(( Jlog2 K + M)KM) for computing the solution on 
M time levels m = 1, . . . , M. 

Notice that the algorithm essentially enlarges possibilities of the corresponding one in 
[25] and also is highly parallelizable. 



4 Numerical experiments 



We have implemented the above algorithm. For numerical experiments, similarly to 
we take p{x^y) = 1, T/q = h = 1 and a simple rectangular potential, i.e. the 
barrier 

(Q for (x, ?/) G (a, 6) X (c, rf) 

v{x,y) = < , g > 0. 

otherwise 

On the other hand, from the numerical point of view, this barrier is not so simple since it is 
discontinuous and consequently the corresponding exact solution is not so smooth. Below 
we choose the fixed (a, b) = (1.6, 1.7) and (c, d) of three different lengths in Examples A, 
B and C. 

Let the initial function be the Gaussian wave package 



^^{x,y) = ^G{x,y) 



exp 



\ iV2k{x 



X 



(0) 



) 



(x 



X 



(0)\2 



)2 + (^_^(0))2 



on 



We take the parameters k = 30 and a = 1/120. 

We solve the initial boundary value problem in the infinite strip M x (0, y), choose the 
computational domain Qx x [0, T] such that (a, b) x (c, d) C Qx and ipc is small enough 
outside Clx' Namely, below X = 3 and Y = 2.8 together with {x^^\ y^^"^) = (1, ^) as well 
as T = = 0.027. We use the uniform meshes in x, y and t with steps respectively 
h 



^ and T 



M' 



We accordingly modify the splitting in potential scheme (2.17)-(2.22) replacing cj/^Uxj 
by cJ/j in (2.17) and (2.19), the boundary conditions (2.20) by 

*'"lw = 0, n=o,K = 0, n=o,K = on uJh 



together with the left discrete TBC at j = similar to the right one (2.21) for S = 5ref. 
On Figure [T] the modulus and the real part of ipc are shown on the computational domain. 
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The normalized barrier (with Q = 1) is situated there as weU, for (c, = (0.7,2.1) (see 
Example B below). 












II 





Figure 1: The modulus (left) and the real part (right) of the initial function together with the normalized 
barrier from Example B 



Example A. We first take (c, rf) = (0, y) and Q = 1500 (previously the same example 
was in fact treated in by a method from [25j). In this case, the wave package is 
divided into two similar refiected and transmitted parts moving in opposite directions 
with respect to the barrier. A little bit surprisingly, the solution is almost the same as 
in the case (c, rf) = (^, ^) in Example B (that is why the corresponding graphs of the 
solution are given below). Namely, for the fine mesh with (J, X, M) = (9600, 512, 4800), 
norms of differences between the pseudo-exact solution for (c, rf) = (0, Y) computed by 
the Crank-Nicolson scheme and one for (c, rf) = ( j, ^) computed by the splitting in 
potential scheme are 

Ec ^ 1.81 • 10"^ El2 ^ 5.02 . 10"^ 

i.e., they are actually small. In this section, we exploit the splitting method with the 
simplest choice V = unless the contrary is explicitly stated. Hereafter Ec and £"7,2 
denote differences/errors in the mesh norms that are uniform in time as well as C (i.e. 
uniform) and in space. 

Notice also that the norms of differences between the pseudo-exact solutions on the 
fine mesh and on the mesh with (J, M) = (1200, 64, 600) are 

£c^2.70•10-^ £;^2 1.65 • lo-^ 

for the Crank-Nicolson scheme and 

2.54 . lO"^ £;^2 1.60- lO"^ 

for the splitting scheme (see also Figure [2] below for more detail), i.e., they are small 
enough and close to each other. 

Example B. Next we consider the barrier with (c, rf) = (j, ^) = (0.7,2.1) (that is 
one-half in length of the first one), for three values of the barrier height in order to get 
qualitatively varying behavior of solutions (compare with [24j). 

First, once again we take the barrier height Q = 1500. We exploit the mesh with 
{J,K,M) = (1200,64,600) so that h = 2.5 • 10-^ 5 = 4.375 • 10"^ and r = 4.5 • 10"^ 
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J 


Er 




Rr 


E 


T 2 




Rr2 


-^time 


300 


0.22 






0.12 










600 


5.62 • 10- 


-2 


3.99 


3.10- 


10- 


-2 


3.94 


1.33 


1200 


1.42 • 10- 


-2 


3.95 


7.91 - 


10- 


-3 


3.92 


1.47 


2 400 


3.77-10- 


-3 


3.77 


2.16- 


10- 


-3 


3.66 


1.63 


4 800 


1.20-10- 


-3 


3.15 


7.65- 


10- 


-4 


2.83 


1.81 



Table 1: Example B for Q = 1500. Errors, ratios of errors and ratios of runtimes in dependence with J, 
for K = 256 and M = 2400 

Though 19 is large, we qualify that such choice is reasonable. In particular, for 

comparison we exploit the above mentioned pseudo-exact solution on the fine mesh with 
(J, X, M) = (9600,512,4800) that all three are 8 times larger and imply the steps h = 
3.125 • 10"^, S ^ 5.469 • 10"^ and r = 5.625 • 10"^. Figure [2] demonsrates the behavior of 
the absolute and relative errors in C and space mesh norms in dependence with time. 



0.03 
0.025 

0.02 
0.015 

0.01 
0.005 



0.009 0.018 0.027 " 0.009 0.018 0.027 

Figure 2: Example B for Q = 1500. The absolute and relative errors in C and norms in dependence with 
time for the numerical solution for (J, M) = (1200, 64, 600) 

The modulus and the real part of the numerical solution are given on Figure [3} for 
the time moments = tut^ m = 180, 300, 420 and 600. In this case, the wave package 
is divided into two rather similar refiected and transmitted parts moving in opposite 
directions with respect to the barrier. 

We continue to study the error behavior in more detail in Tables [T} [2] and (S) They 
contain errors of the solutions to the splitting method for increasing J, K and M respec- 
tively (for sufficiently large values of two other numbers). The associated ratios Rc and 
Rl2 of the sequential errors are also put there. They are rather close to 4 excluding the 
last rows that means almost the second order of convergence with respect to each of J, K 
and M, both in C and mesh norms. The deterioration of Rc and Rl2 in the last rows 
is explained by more essential inffuence of the errors due to the chosen discretization in 
other directions. 

The last columns in the tables contain also the respective ratios of runtimes. One can 
see that they all are close to 2 when any of J, X or M increases twice. 

In Table [i] we put C and errors for some selected values of J, X and M. They 
all decrease monotonically as J, X or M increase. We also compare there the numer- 



Absolute error Relative error 
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K 


Ec 




Rc 


El- 




Rl- 


-^time 


16 


0.15 






6.62 • 10- 


-2 






32 


3.43 • 10" 


-2 


4.36 


2.34-10- 


-2 


2.84 


2.13 


64 


8.65 • 10" 


-3 


3.96 


6.58 - 10- 


-3 


3.55 


1.91 


128 


2.29 • 10- 


-3 


3.79 


1.81 - 10- 


-3 


3.63 


1.86 


256 


1.20-10- 


-3 


1.91 


7.65 - 10- 


-4 


2.37 


1.87 



Table 2: Example B for Q = 1500. Errors, ratios of errors and ratios of runtimes in dependence with 
for J = 4800 and M = 2400 



M 


Ec 




Rc 


El- 




Rl- 


-^time 


150 


0.17 






9.10-10" 


-2 






300 


4.39 - 10- 


-2 


3.91 


2.37-10 


-2 


3.84 


2.1 


600 


1.12-10- 


-2 


3.9 


6.20 - 10 


-3 


3.83 


2.02 


1200 


3.16-10- 


-3 


3.56 


1.83-10" 


-3 


3.39 


2.01 


2 400 


1.20-10- 


-3 


2.64 


7.65 - 10 


-4 


2.39 


2.09 



Table 3: Example B for Q = 1500. Errors, ratios of errors and ratios of runtimes in dependence with M, 
for J = 4800 and K = 256 

ical solutions of the splitting method with 1/ = and V{x) = Qx{^)^ where x(^) is 
the characteristic function of the interval (a, 6); two last columns of the table contain 
percentages 

( ^^}y=o _ i\ . 100%, P^. := (44^ - • 100%. 

One can see that the second choice V = Qx also works but the first one V = mostly 
leads to better results. 

In addition, for the fine mesh with {J,K,M) =^9600,512,4800) the norms of differ- 
ences between the solutions for these two different V are 

Ec ^ 3.32 • 10-^ ^ 1.04 • 10-^ 

i.e., they are very small. 



J 


K 


M 


Ec 




El- 




Pc 


Pl- 


1200 


64 


600 


2.54-10- 


-2 


1.60-10- 


-2 


-5.65 


-3.08 


1200 


128 


600 


2.42 • 10- 


-2 


1.37-10- 


-2 


-6.6 


-3.78 


1200 


64 


1200 


1.83 • 10- 


-2 


1.21 - 10- 


-2 


0.36 


-0.91 


1200 


128 


1200 


1.63 • 10- 


-2 


9.40 • 10- 


-3 


-2.49 


-1.33 


2 400 


64 


600 


1.64-10- 


-2 


1.09 - 10- 


-2 


-5.48 


-3.99 


2 400 


128 


600 


1.38 • 10- 


-2 


8.05 - 10- 


-3 


-11.39 


-6.27 


2 400 


64 


1200 


1.03 • 10- 


-2 


7.73 • 10- 


-3 


1.01 


-0.96 


2 400 


128 


1200 


6.00 • 10- 


-3 


3.82 - 10- 


-3 


-6.3 


-3.02 



Table 4: Example B for Q = 1500. Errors of the numerical solutions for y = and percentages of their 
changes when taking V = Qx 
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J K M Ec Ec, rel ^L^rel 

2 400 64 600 8.61 • 10"^ 5.29 • 10"^ 2.48 • 10-^ 2.90 • 10-^ 

1200 128 600 1.38-10-2 6.37-10-^ 3.56 • 10"^ 2.82 • 10"^ 

1200 64 1200 1.22-10-2 5.30-10-^ 2.82-10-2 2.33 - 10-^ 



Table 5: Example B for Q = 4000. The change in numerical solution when J ^ K oi M increases twice 



J K M Ec El2 Ec, rel ^L^rel 

2 400 512 600 1.06 - 10-^ 5.26 - 10-^ 3.16 - 10-^ 2.32 - 10-^ 

1200 1024 600 9.64-10-^ 6.70-10-^ 3.92 - 10-^ 4.38 - 10-^ 

1200 512 1200 8.22-10-^ 4.67-10-^ 2.66 - 10-^ 2.06 - 10-^ 



Table 6: Example C. The change in numerical solution when J ^ K oi M increases twice 

Second, we take a less barrier height Q = 1000. This situation is simpler from the 
numerical point of view. The numerical results are demonstrated on Figure [4] for the 
same time moments and mesh. Now the wave package goes through the barrier with an 
essentially less reflection. 

Third, let Q = 4000 be rather large. On Figure [5] the numerical solution is represented 
for the same time moments and mesh. Here the main part of the wave is reflected from 
the barrier and then moves in the opposite direction along the x axis. 

To check the approximate solution in this case, we compute how the numerical solution 
changes when any of J, X or M increases twice, see Table [5} where the corresponding 
absolute and relative errors in C and L2 norms are given. The relative errors Ec^rei and 
EL'^^rei are deflned as the maximal in time relative C and L2 mesh errors in space (in joint 
nodes). One can see that all the errors are small enough. 

We emphasize that on all Figures |3| |4] and [5| the last two graphs exhibit complete 
absence of the spurious reflections from the artiflcial left and right boundaries due to 
exploiting of the discrete TBCs there. 

Example C. We also treat the case of a very short barrier with (c, rf) = — |^ + 
^) = (1.3125,1.4875) and once again having the height Q = 1500; this barrier looks like 
a column. The numerical solution is represented on Figure [6| for the time moments 
tjn = mr^ m = 180, 240, 300 and 360, together with the normalized barrier. We use the 
mesh with (J, X, M) = (1200, 512, 600), i.e. for the same J and M as on the above flgures 
but for notably larger K. Now in contrast to the previous examples, the transmitted part 
of the wave package is separated into two pieces. 

To check the approximate solution in this case, we compute how the numerical solution 
changes when J, K or M increases twice, see Table |6| where the corresponding absolute 
and relative errors are given. One can see that all of them are small enough once again. 

We call attention to the essentially more complicated behavior of the real part of the 
solution compared to its modulus in all Examples A-C. Comparing C and errors, one 
can see that though C errors are mainly larger, their behavior is rather similar that is not 
so obvious a priori taking into account the oscillatory type of the solutions in space and 
time. 

In general, the above practical error analysis indicates the good error properties of the 
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splitting in potential scheme. 

Finally, note that clearly both the rectangular form of the barrier and the specific 
choice of the initial function are inessential to apply efficiently the splitting method. 
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Figure 3: Example B for Q — 1500. The modulus and the real part of the numerical solution m — 

180, 300, 420 and 600 for (J, K, M) = (1200, 64, 600) 
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Figure 4: Example B for Q — 1000. The modulus and the real part of the numerical solution m = 

180,300,420 and 600 for {J,K,M) = (1200,64,600) 
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Figure 5: Example B for Q = 4000. The modulus and the real part of the numerical solution m 
180, 300, 420 and 600 for (J, M) = (1200, 64, 600) 
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Figure 6: Example C. The modulus and the real part of the numerical solution m 
for {J,K,M) = (1200,512,600) 



180,240,300 and 360 
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